Prep

setwd('/project/CRI/DeBerardinis_lab/lcai/NSCLC_Metabolism/scripts/final/')
.libPaths(c(.libPaths(),'~/R/x86_64-unknown-linux-gnu-library/3.2/'))
# libraries
library(GGally)
library(openxlsx)
library(Hmisc)
library(corrplot)
library(RColorBrewer)
library(GSVA)
library(ppcor)
library(reshape2)
library(gtable)
library(grid)
library(gridExtra)
library(factoextra)
library(ggrepel)
library(heatmap3)
library(stringr)
library(plotrix)
library(mclust)
library(MASS)
library(sjPlot)
library(ComplexHeatmap)
library(ClassComparison)
library(pracma)
library(ggpubr)
library(cowplot)

# load data 
c2cgp<-readRDS('data/c2cgp.rds')
c2cp<-readRDS('data/c2cp.rds')
dat<-readRDS('data/metabolic_profiling_data.rds')
ge<-readRDS('data/illumina_gene_level.rds')
common<-intersect(rownames(ge),rownames(dat))
ge<-ge[common,]
# ssGSEA.results<-gsva(expr = t(ge),gset.idx.list = c2cgp,method="ssgsea",rnaseq=F)
# The line above is commented because it takes too long to run, the pre-calculated results were saved and loaded below
ssGSEA.results<-readRDS('data/ssGSEA_result.rds')
sig.ne <- read.xlsx('data/Zhang_2018_NE_signature.xlsx')
pc<-readRDS('data/34glucose_PC.rds')
mutations<-readRDS('data/mutations.rds')
dat.invivo<-readRDS('data/invivo_data.rds')
mut<-readRDS('data/invivo_mutation.rds')
lkb1<-readRDS('data/LKB1_OE.rds')
S4.data<-readRDS('data/Table_S4.rds')
AUC<-read.xlsx('data/table6_chemical_screen_data.xlsx',sheet = 'AUC')
emt.signature<-read.csv('data/EMT_signature.csv')
emt.signature<-unique(emt.signature$Gene.Symbol)
emt.signature<-gsub(' ','',fixed = T,emt.signature)
pem<-readRDS('data/pem_validation.rds')
wb.misty<-loadWorkbook('data/Misty_xenograft.xlsx')
wb<-loadWorkbook('data/All Metabolic data - 2.16.2015R1.xlsx')
CDH1.data<-readRDS('data/CDH1_multiomics.rds')
load('data/additional_serine_data.RData')
col2<-colorRampPalette(rev(brewer.pal(9,'RdBu')))(100)

# Shared functions
cor_fun <- function(data, mapping, method="pearson", ndp=2, sz=5, stars=TRUE, ...){
  xVal<-mapping_string(mapping$x)
  yVal<-mapping_string(mapping$y)
  data <- na.omit(data[,c(xVal,yVal)])

  x<-data[,xVal]
  y<-data[,yVal]
  
  corr <- cor.test(x, y, method=method)
  est <- corr$estimate
  
  if(stars){
    stars <- c("***", "**", "*", "")[findInterval(corr$p.value, c(0, 0.001, 0.01, 0.05, 1))]
    lbl <- paste0(round(est, ndp), stars)
  }else{
    lbl <- round(est, ndp)
  }
  
  ggplot(data=data, mapping=mapping) + 
    annotate("text", x=mean(x), y=mean(y), label=lbl)+
    theme(panel.grid = element_blank())
}

col.grad<-function(x){
  scales::alpha(colorRampPalette(rev(brewer.pal(5,'RdBu')))(length(x)),alpha = 0.5)[rank(x)]
}
colorlegend <- function(
  colbar,
  labels,
  at = NULL,
  xlim = c(0, 1),
  ylim = c(0, 1),
  vertical = TRUE,
  ratio.colbar = 0.4,
  lim.segment = "auto", # NOTE: NULL treated as "auto"
  align = c("c", "l", "r"),
  addlabels = TRUE,
  ...)
{
  if (is.null(at) && addlabels) {
    at <- seq(0L, 1L, length = length(labels))
  }
  
  if (is.null(lim.segment) || lim.segment == "auto") {
    lim.segment <- ratio.colbar + c(0, ratio.colbar * .2)
  }
  
  if (any(at < 0L) || any(at > 1L)) {
    stop("at should be between 0 and 1")
  }
  
  if (length(lim.segment) != 2) {
    stop("lim.segment should be a vector of length 2")
  }
  
  if (any(lim.segment < 0L) || any(lim.segment > 1L)) {
    stop("lim.segment should be between 0 and 1")
  }
  
  align <- match.arg(align)
  xgap <- diff(xlim)
  ygap <- diff(ylim)
  len <- length(colbar)
  rat1 <- ratio.colbar
  rat2 <- lim.segment
  
  if (vertical) {
    
    at <- at * ygap + ylim[1]
    yyy <- seq(ylim[1], ylim[2], length = len + 1)
    
    rect(rep(xlim[1], len), yyy[1:len],
         rep(xlim[1] + xgap * rat1, len), yyy[-1],
         col = colbar, border = colbar)
    rect(xlim[1], ylim[1], xlim[1] + xgap * rat1, ylim[2], border = "black")
    segments(xlim[1] + xgap * rat2[1], at, xlim[1] + xgap * rat2[2], at)
    
    if (addlabels) {
      pos.xlabel <- rep(xlim[1] + xgap * max(rat2, rat1), length(at))
      switch(align,
             l = text(pos.xlabel, y = at, labels = labels, pos = 4, ...),
             r = text(xlim[2],    y = at, labels = labels, pos = 2, ...),
             c = text((pos.xlabel + xlim[2]) / 2, y = at, labels = labels, ...),
             stop("programming error - should not have reached this line!")
      )
    }
  } else {
    
    at <- at * xgap + xlim[1]
    xxx <- seq(xlim[1], xlim[2], length = len + 1)
    
    rect(xxx[1:len], rep(ylim[2] - rat1 * ygap, len),
         xxx[-1], rep(ylim[2], len),
         col = colbar, border = colbar)
    rect(xlim[1], ylim[2] - rat1 * ygap, xlim[2], ylim[2], border = "black")
    segments(at, ylim[2] - ygap * rat2[1], at, ylim[2] - ygap * rat2[2])
    
    if (addlabels) {
      pos.ylabel <- rep(ylim[2] - ygap * max(rat2, rat1), length(at))
      switch(align,
             l = text(x = at, y = pos.ylabel, labels = labels, pos = 1, ...),
             r = text(x = at, y = ylim[1],    labels = labels, pos = 2, ...),
             c = text(x = at, y = (pos.ylabel + ylim[1]) / 2, labels = labels, ...),
             stop("programming error - should not have reached this line!")
      )
    }
  }
}

addLegend<-function(x,name){
  colbar<-col.grad(sort(x))
  labels<-cl.lim<-round(range(x,na.rm = T),1);cl.length<-100
  cl.align.text = "c"
  par(mar=c(0.5,3.2,2,3.2),xpd=T)
  plot.new()
  title(name,line = 0)
  colorlegend(colbar = colbar, labels = labels,
              offset = 1, ratio.colbar = 0.5, cex = 1,
              vertical = F,
              align = cl.align.text)
}

tri_plot<-function(x){
  addLegend(x[,3],name=names(x)[3])
  par(mar=c(3,3,0.2,3))
  plot(x[,1],x[,2],col=col.grad(x[,3]),pch=19,xlab=names(x)[1],ylab=names(x)[2],mgp=c(2,1,0),oma=c(2,0,2,0))
  
  addLegend(x[,1],name=names(x)[1])
  par(mar=c(3,3,0.2,3))
  plot(x[,2],x[,3],col=col.grad(x[,1]),pch=19,xlab=names(x)[2],ylab=names(x)[3],mgp=c(2,1,0),oma=c(2,0,2,0))
  
  addLegend(x[,2],name=names(x)[2])
  par(mar=c(3,3,0.2,3))
  plot(x[,3],x[,1],col=col.grad(x[,2]),pch=19,xlab=names(x)[3],ylab=names(x)[1],mgp=c(2,1,0),oma=c(2,0,2,0))
}

Isotopologue.Distribution<-function(Key,Name){
  dat.key<-dat[,grep(Key,grep("m[0-9]$",names(dat),value = T))]
  dat.key<-setNames(melt(dat.key[complete.cases(dat.key),]),c('Isotopologue','% Enrichment'))
  dat.key$tracer<-factor(ifelse(grepl('Q',dat.key$Isotopologue),'Glutamine','Glucose'),levels=c('Glucose','Glutamine'))
  dat.key$timepoint<-factor(ifelse(grepl('24',dat.key$Isotopologue),'24h','6h'),levels=c('6h','24h'))
  dat.key$Isotopologue<-sapply(dat.key$Isotopologue,FUN = function(x){
    tmp<-unlist(strsplit(split = 'm',as.character(x)))
    return(paste0('m+',tmp[length(tmp)]))
  })
  ggplot(dat.key, aes(x=Isotopologue, y=`% Enrichment`)) + 
    geom_violin(trim = T,fill='#6e016b',alpha=0.5)+
    geom_boxplot(width=0.1,outlier.shape = NA)+
    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
    theme_bw()+
    facet_grid(tracer~timepoint)+
    ggtitle(paste0(Name,' Labeling'))
}

Fig 1

1B

ggpairs(dat[,c('Glc','Lac','Gln','Glu')],
           upper=list(continuous=cor_fun),
           lower=list(continuous = wrap("points", alpha = 0.5,colour='#2c7fb8')))+
  theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())

1C-F

LacGlc_ssGSEA<-merge(dat[,'Lac/Glc',drop=F],t(ssGSEA.results),by = 0,all = F)
LacGlc_ssGSEA.r<-apply(LacGlc_ssGSEA[,-c(1:2)],2,FUN=function(x) cor.test(x,LacGlc_ssGSEA$`Lac/Glc`)$estimate)
lacGlc_hypoxia<-LacGlc_ssGSEA.r[grep('HYPOXIA',names(LacGlc_ssGSEA.r))]
lacGlc_hypoxia<-lacGlc_hypoxia[grep('DN',names(lacGlc_hypoxia),invert = T)]

layout(matrix(c(1,2,4,3,2,4),byrow = T,ncol = 3,nrow = 2),widths=c(3,2,2), heights=c(1,1))

plot(density(lacGlc_hypoxia,cut = range(lacGlc_hypoxia),bw = 0.05),col='#810f7c',lwd=2,
     main='Correlation between Lac/Glc and Hypoxia Gene Signature Enrichment Scores',
     cex.main=0.9)
points(lacGlc_hypoxia,jitter(rep(0.55,length(lacGlc_hypoxia)),amount = 0.1),col=scales::alpha('#810f7c',0.4),pch=19)

tmp<-LacGlc_ssGSEA[,c('WINTER_HYPOXIA_UP','Lac/Glc')]
tmp.cor<-rcorr(as.matrix(tmp))
plot(tmp,
     main=paste0('r = ',round(LacGlc_ssGSEA.r['WINTER_HYPOXIA_UP'],2),', pv = ',signif(tmp.cor$P[1,2],2)),
     pch=19,col=scales::alpha(colour = '#d7301f',alpha = 0.5),
     xlab='WINTER_HYPOXIA_UP enrichment score')

winter_hypoxia_r<-rcorr(cbind(dat[common,]$`Lac/Glc`,ge[common,intersect(c2cgp$WINTER_HYPOXIA_UP,colnames(ge))]))$r[1,-1]
plot(density(winter_hypoxia_r,bw = 0.05,cut = range(winter_hypoxia_r)),col='#810f7c',lwd=2,
     main='Correlation between Lac/Glc and WINTER_HYPOXIA_UP genes',cex.main=0.9)
set3<-winter_hypoxia_r[!names(winter_hypoxia_r)%in%c(c2cp$REACTOME_GLYCOLYSIS,'LDHA')]
points(set3,jitter(rep(0.55,length(set3)),amount = 0.1),col=scales::alpha('#810f7c',0.1),pch=19)
set1<-winter_hypoxia_r[names(winter_hypoxia_r)%in%c2cp$REACTOME_GLYCOLYSIS]
points(set1,jitter(rep(0.55,length(set1)),amount = 0.1),col=scales::alpha('black',0.8),pch=19)
set2<-winter_hypoxia_r[names(winter_hypoxia_r)=='LDHA']
points(set2,jitter(rep(0.55,length(set2)),amount = 0.1),col=scales::alpha('red',0.8),pch=19)
legend('topright',legend = c('REACTOME_GLYCOLYSIS','LDHA','other'),bty='n',bg = 'transparent',cex = 0.8,pch=19,
       col=scales::alpha(alpha = c(0.8,0.8,0.1),colour = c('black','red','#810f7c')))

comp.score <- function(dat, sig) {
  ind <- match(sig$Symbol,rownames(dat))
  ind <- ind[!is.na(ind)]
  sig <- sig[sig$Symbol %in% rownames(dat), ]
  score1 <- apply(dat[ind, ], 2, function(d) cor.test(sig[, 'NE.class.mean.expression'], d)$estimate)
  score2 <- apply(dat[ind, ], 2, function(d) cor.test(sig[, 'Non-NE.class.mean.expression'], d)$estimate)
  score <- (score1 - score2)/2
  score
}

ne_score<-comp.score(t(ge),sig = sig.ne)
tmp<-data.frame(ne_score,dat[common,'Lac/Glc'])
tmp<-tmp[complete.cases(tmp),]
tmp.cor<-rcorr(as.matrix(tmp))
plot(tmp,main=paste0('r = ',round(tmp.cor$r[1,2],2),', pv = ',signif(tmp.cor$P[1,2],2)),
     pch=19,col=scales::alpha(colour = '#d7301f',alpha = 0.5),
     xlab='neuroendocrine score',ylab='Lac/Glc')

1G

corrplot(rcorr(as.matrix(data.frame('Glu/Gln'=dat[common,c('Glu/Gln')],ge[,c('GLS','GLS2')],check.names = F)))$r,cl.pos = 'n',
         col=col2,addCoef.col = "black",diag = F,type = 'lower',method = 'color',tl.col='black')

1H

corrplot(rcorr(as.matrix(dat[,c('Glc','Lac','Lac/Glc','Gln','Glu','Glu/Gln',
                                'Day3/Day1','Day5/Day1',
                                'Day1-G','Day3-G','Day5-G','Day1-Q','Day3-Q','Day5-Q')]))$r,
         col=col2,diag = F,type = 'upper',method = 'color',addCoef.col = "black",tl.col='black')

1I top

make_cormat<-function(x){
  corMat<-cor(x)
  pcorMat<-pcor(x)$estimate
  colnames(pcorMat)<-rownames(pcorMat)<-colnames(corMat)
  return(list(corMat,pcorMat))
}

plot_cor<-function(corMat,pcorMat,cex=1){
  par(mfrow=c(1,2))
  corrplot(corMat,col=col2,addCoef.col = "black",method = 'color',cl.pos = 'n',
           type = 'upper',tl.cex = cex,tl.col = 'black',diag = F,title = 'pairwise correlation',mar=c(2,2,4,2))
  corrplot(pcorMat,col=col2,addCoef.col = "black",method = 'color',cl.pos = 'n',
           type = 'upper',tl.cex = cex,tl.col = 'black',diag = F,title = 'partial correlation',mar=c(2,2,4,2))
}

par(mfrow=c(1,2))
tmp<-make_cormat(na.omit(dat[,c('Lac','Glc','Day5-G')]));plot_cor(tmp[[1]],tmp[[2]])

1J top

par(mfrow=c(1,2))
tmp<-make_cormat(na.omit(dat[,c('Lac','Glc','Gln')]));plot_cor(tmp[[1]],tmp[[2]])

1I-J bottom

layout(matrix(1:12,ncol = 6,byrow = F),widths=rep(2,6), heights=c(2, 5))
tri_plot(dat[,c('Glc','Lac','Day5-G')])
tri_plot(dat[,c('Glc','Lac','Gln')])

Fig 2

2C

Isotopologue.Distribution(Key = 'Cit',Name='Citrate')

2D-E

c1<-rcorr(as.matrix(dat[,c('CitG6m0','CitQ6m5')]))
c2<-rcorr(as.matrix(dat[,c('CitG6m2','CitQ6m4')]))
marrangeGrob(list(ggplot(dat[,c('CitG6m0','CitQ6m5')],mapping = aes(x=CitG6m0,y=CitQ6m5))+
                    geom_point(colour='#6e016b',alpha=0.5,size=2.5)+theme_bw()+
                    ggtitle(paste0("Glutamine dependent reductive carboxylation\nr = ",round(c1$r[1,2],2),', pv = ',signif(c1$P[1,2],2))),
                  ggplot(dat[,c('CitG6m2','CitQ6m4')],mapping = aes(x=CitG6m2,y=CitQ6m4))+
                    geom_point(colour='#6e016b',alpha=0.5,size=2.5)+theme_bw()+
                    ggtitle(paste0("Glutamine dependent anaplerosis\nr = ",round(c2$r[1,2],2),', pv = ',signif(c2$P[1,2],2)))), 
             nrow=2, ncol=1, top=NULL)

2F

metabolite<-sapply(colnames(dat),FUN = function(x) unlist(strsplit(x,split = "m[0-9]"))[1])
select<-names(table(metabolite))[table(metabolite)>1]
dat.sub<-dat[,metabolite%in%select]
metabolite.sub<-sapply(colnames(dat.sub),FUN = function(x) unlist(strsplit(x,split = "m[0-9]"))[1])
G6<-grep('G6',unique(metabolite.sub),value = T)
Q6<-grep('Q6',unique(metabolite.sub),value = T)
G24<-grep('G24',unique(metabolite.sub),value = T)
Q24<-grep('Q24',unique(metabolite.sub),value = T)
nutrient.contribution<-lapply(list(G6,Q6,G24,Q24),FUN=function(c){
  setNames(as.data.frame(do.call(cbind,lapply(c,FUN=function(x){
    apply(dat.sub[,metabolite.sub==x],1,FUN=function(y){
      sum((0:(length(y)-1))*(y/100))/(length(y)-1)
    })
  }))),nm = c)
})
names(nutrient.contribution)<-c('G6','Q6','G24','Q24')
tmp6<-setNames(data.frame(melt(nutrient.contribution$G6[,gsub('Q','G',colnames(nutrient.contribution$Q6))]),melt(nutrient.contribution$Q6)[,2]),
               nm = c('experiment','[U-13C]glucose','[U-13C]glutamine'))
tmp24<-setNames(data.frame(melt(nutrient.contribution$G24[,gsub('Q','G',colnames(nutrient.contribution$Q24))]),melt(nutrient.contribution$Q24)[,2]),
                nm = c('experiment','[U-13C]glucose','[U-13C]glutamine'))
tmp<-data.frame(time=rep(c('6h','24h'),c(nrow(tmp6),nrow(tmp24))),rbind(tmp6,tmp24),check.names = F)
tmp$experiment<-gsub('G',' ',tmp$experiment)
tmp$experiment<-gsub("$","h",tmp$experiment)
tmp$experiment<-factor(tmp$experiment)
tmp$experiment<-factor(tmp$experiment,levels=levels(tmp$experiment)[c(2,4,8,6,1,3,7,5)])

tmp$sum<-tmp$`[U-13C]glucose`+tmp$`[U-13C]glutamine`
ggplot(tmp,aes(x=sum))+geom_density()+facet_wrap(~experiment,nrow = 2)+theme_bw()+xlab("sum of glucose and glutamine fractional carbon contribution")

2G

h<-hclust(as.dist(1-abs(rcorr(as.matrix(dat),type = 'pearson')$r)),method = 'ward.D2')
fviz_dend(h, k = 23, k_colors = "jco",main = 'Clustering Metabolic Features',
          type = "phylogenic", repel = TRUE)

Fig 3

3C

colnames(pc)[1]<-'cell line'
pc$select<-rep(NA,nrow(pc))
pc$select[pc$`cell line`%in%c('H920','PC-9','H2444')]<-'low'
pc$select[pc$`cell line`%in%c('HCC515','H1792','H1648')]<-'high'
ct<-cor.test(pc$`m+1 Citrate [3,4-13C]`,pc$`m+3 malate [U-13C]`)
ggplot(pc,aes(x=`m+1 Citrate [3,4-13C]`,y=`m+3 malate [U-13C]`))+
  geom_point(aes(colour=select),alpha=0.7,size=2.5)+ 
  geom_text_repel(data=subset(pc,!is.na(select)),aes(label=`cell line`),size = 3.5)+
  scale_colour_discrete(name="Cell Lines",
                        breaks=c("high", "low", NA),
                        labels=c("high", "low", "unselected"))+
  ggtitle(paste("pyruvate carboxylase dependent anaplerosis\nr = ",
                round(ct$estimate,2),", pv = ",signif(ct$p.value,2),sep = ''))+
  
  labs(x=expression('m+1 citrate from [3,4-'^{13}*'C] glucose labeling'),y=expression('m+3 malate from [U-'^{13}*'C] glucose labeling'))+
  theme_light()

3D was made by Pei-Husan Chen with Prism GraphPad (see data/Figure_3D.pzfx)

Fig 4

4A

# KRAS mutations
KRAS<-mutations[,'KRAS']
KRAS.mut<-which(KRAS!="")
KRAS.wt<-which(KRAS=="")
KRAS.ms<-grep('MS',KRAS,ignore.case = F)
aa_change<-str_extract(pattern = "[pP].[ ]?[A-Z][0-9]{1,5}[A-Z]", string = KRAS[KRAS.ms])
aa_site<-as.numeric(str_extract(pattern = "[0-9]+",string = aa_change))
KRAS.ms12<-KRAS.ms[aa_site==12]
KRAS.ms13<-KRAS.ms[aa_site==13]
KRAS.ms61<-KRAS.ms[aa_site==61]
# EGFR mutations
EGFR<-mutations[,'EGFR']
EGFR.mut<-which(EGFR!="")
exon19.del<-str_extract(pattern = "[A-Z]7[0-9]+[-_][A-Z]7[0-9]+[ ]*del",string = EGFR)
EGFR.wt<-which(EGFR=='')
EGFR.exon19del<-which(!is.na(exon19.del))
# STK11 mutations
STK11<-mutations[,'STK11']
STK11.mut<-which(STK11!="")
STK11.wt<-which(STK11=="")
STK11.ms<-grep('MS',STK11,ignore.case = F)
STK11.ns<-grep('NS',STK11,ignore.case = F)
STK11.ss<-grep('SS',STK11,ignore.case = F)
STK11.fs<-grep('fs;',STK11,ignore.case = F)
STK11.ms<-setdiff(STK11.ms,STK11.fs)
KRAS_STK11<-intersect(KRAS.mut,STK11.mut)
select.mut<-do.call(cbind,lapply(list(KRAS.wt,KRAS.mut,KRAS.ms12,KRAS.ms13,KRAS.ms61,EGFR.wt,EGFR.mut,EGFR.exon19del,STK11.wt,STK11.mut,STK11.ms,STK11.ns,STK11.fs,STK11.ss), FUN=function(x) (1:nrow(mutations))%in%x))
colnames(select.mut)<-c('KRAS.wt','KRAS.mut','KRAS.ms12','KRAS.ms13','KRAS.ms61','EGFR.wt','EGFR.mut','EGFR.exon19del','STK11.wt','STK11.mut','STK11.ms','STK11.ns','STK11.fs','STK11.ss')
specific.mut.mat<-matrix(NA,ncol=ncol(dat),nrow = ncol(select.mut)-3,dimnames = list(grep('wt',colnames(select.mut),invert = T,value = T),colnames(dat)))
for(i in 1:nrow(specific.mut.mat)){
  pick<-grep('wt',colnames(select.mut),invert = T,value = T)[i]
  pick.gene<-unlist(strsplit(pick,split = '.',fixed = T))[1]
  pick.ctrl<-paste0(pick.gene,'.wt')
  for(j in 1:ncol(dat)){
    specific.mut.mat[i,j]<-t.test(dat[select.mut[,pick.ctrl],j],dat[select.mut[,pick],j])$p.value
  }
}
library(RColorBrewer)
KRAS<-ifelse(select.mut[,'KRAS.wt'],'white','brown')
KRAS.col<-brewer.pal(3,'Set1')
KRAS[select.mut[,'KRAS.ms12']]<-KRAS.col[1]
KRAS[select.mut[,'KRAS.ms13']]<-KRAS.col[2]
KRAS[select.mut[,'KRAS.ms61']]<-KRAS.col[3]

EGFR<-ifelse(select.mut[,'EGFR.wt'],'white','brown')
EGFR.col<-brewer.pal(3,'Set2')
EGFR[select.mut[,'EGFR.exon19del']]<-EGFR.col[1]

STK11<-ifelse(select.mut[,'STK11.wt'],'white','brown')
STK11.col<-brewer.pal(4,'Set3')
STK11[select.mut[,'STK11.ms']]<-STK11.col[1]
STK11[select.mut[,'STK11.ns']]<-STK11.col[2]
STK11[select.mut[,'STK11.fs']]<-STK11.col[3]
STK11[select.mut[,'STK11.ss']]<-STK11.col[4]

KRAS_col<-ifelse(select.mut[,'KRAS.wt'],'#f1eef6','#dd1c77')
KRAS.ms12_col<-ifelse(!select.mut[,'KRAS.ms12'],'#f1eef6',"#df65b0")
KRAS.ms13_col<-ifelse(!select.mut[,'KRAS.ms13'],'#f1eef6',"#df65b0")
KRAS.ms61_col<-ifelse(!select.mut[,'KRAS.ms61'],'#f1eef6',"#df65b0")
EGFR_col<-ifelse(select.mut[,'EGFR.wt'],'#f1eef6','#dd1c77')
EGFR.exon19del_col<-ifelse(!select.mut[,'EGFR.exon19del'],'#f1eef6','#df65b0')
STK11_col<-ifelse(select.mut[,'STK11.wt'],'#f1eef6','#dd1c77')
KRAS_STK11_col<-ifelse(!(select.mut[,'STK11.mut'] & select.mut[,'KRAS.mut']),'#f1eef6','#980043')
col.mat<-cbind(EGFR_col,EGFR.exon19del_col,KRAS_col,KRAS.ms12_col,KRAS.ms13_col,KRAS.ms61_col,STK11_col,KRAS_STK11_col)
colnames(col.mat)<-gsub('_col','',fixed = T,colnames(col.mat))
heatmap3(dat[complete.cases(dat[,1:7]),1:7],
         dist=dist,Colv = NA,col = colorRampPalette(c('#f1eef6','red'))(1000),method = 'ward.D2',balanceColor = F,scale='none',
         RowSideColors = col.mat,labCol = paste('m+',0:6),labRow = NA)
legend('topleft',lty=1,lwd=3,legend=c('co-mutation','gene-specific mutation','site-specific mutation','others'),
       col=c('#980043','#dd1c77','#df65b0','#f1eef6'),cex = 0.7,bty='n',inset = c(0.1,-0.2),xpd = T)

4B-E

ecdf_mut<-function(mut,col_choice,mut_name,rest){
  k<-ks.test(dat$CitG6m0[setdiff(1:nrow(dat),mut)],dat$CitG6m0[mut])
  plot(ecdf(dat$CitG6m0[mut]), col=c('#df65b0','#980043')[col_choice], main=NA,
       xlim=range(dat$CitG6m0,na.rm = T),xlab='CitG6m0',ylab='Cumulative Probability')
  plot(ecdf(dat$CitG6m0[setdiff(1:nrow(dat),mut)]), col='#ccc6d6', add=T,main=NA)
  legend(34,0.3,pch=19,col=c(c('#df65b0','#980043')[col_choice],'#ccc6d6'),legend = c(mut_name,rest),bty='n')
  legend(8,0.99,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
}
par(mfrow=c(4,1),mar=c(3,3,1,1),mgp=c(2,1,0))
ecdf_mut(mut = EGFR.exon19del,col_choice = 1,mut_name = 'EGFR.exon19del',rest='Others')
ecdf_mut(mut = KRAS.ms13,col_choice = 1,mut_name = 'KRAS.ms13',rest='Others')
ecdf_mut(mut = KRAS.ms61,col_choice = 1,mut_name = 'KRAS.ms61',rest='Others')
ecdf_mut(mut = KRAS_STK11,col_choice = 2,mut_name = 'KRAS_STK11',rest='Others')
x<-dat$CitG6m0[match(c('A549','HCC44','H460'),rownames(dat))]
y<-(1/length(KRAS_STK11))*match(x,sort(dat$CitG6m0[KRAS_STK11]))
text(x+2,y,c('A549','HCC44','H460'),pos = 4,col='#980043')

4F

produce.mut<-function(x,gene){
  y<-x
  y[grepl(gene,x)]<-'Mut'
  y[grepl('ND',x)]<-'WT'
  y[grepl('NT',x)]<-NA
  y[!y%in%c('Mut','WT',NA)]<-'WT'
  y
}
mut$EGFR<-produce.mut(mut$Mutation,'EGFR')

dat.tumor<-dat.invivo[!grepl('adjacent',dat.invivo$Tissue.fragment,ignore.case = T),]
dat.normal<-dat.invivo[grepl('adjacent',dat.invivo$Tissue.fragment,ignore.case = T),]
dat.tumor[,3:9]<-dat.tumor[,3:9]*100
dat.normal[,3:9]<-dat.normal[,3:9]*100

ecdf_mut<-function(mut,col_choice,mut_name,rest,feature){
  k<-ks.test(dat[,feature][setdiff(1:nrow(dat),mut)],dat[,feature][mut])
  plot(ecdf(dat[,feature][mut]), col=c('#df65b0','#980043')[col_choice], main=NA,
       xlim=range(dat[,feature]),xlab=feature,ylab='Cumulative Probability')
  plot(ecdf(dat[,feature][setdiff(1:nrow(dat),mut)]), col='#ccc6d6', add=T,main=NA)
  legend(34,0.3,pch=19,col=c(c('#df65b0','#980043')[col_choice],'#ccc6d6'),legend = c(mut_name,rest),bty='n')
  legend(8,0.98,bty='n',legend = paste0('pv =',signif(k$p.value,2)))
}

par(oma = c(6,6,4,4) + 0.1,
    mar = c(0,0,1,1) + 0.1)
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE), 
       widths=c(5,5,2))
merge.dat<-merge(dat.tumor,mut,by.x='Pt.Code',by.y='X1',all = F)
k<-ks.test(merge.dat$`M+0`[merge.dat$EGFR=='Mut'],merge.dat$`M+0`[merge.dat$EGFR=='WT'])
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='WT']), col='#ccc6d6', main='Tumor')
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']), pch=19, cex=1.5,add=T,
     xlim=c(50,100),main=NA,col="#df65b0",col.01line = "#df65b0")
lines(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']),col="#df65b0",pch=19,cex=0.6,
      col.points=brewer.pal(name = 'Set1',n = 6)[match(merge.dat$Pt.Code[merge.dat$EGFR=='Mut'],
                                                       unique(merge.dat$Pt.Code[merge.dat$EGFR=='Mut']))])

legend(50,1,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
merge.dat<-merge(dat.normal,mut,by.x='Pt.Code',by.y='X1',all = F)
k<-ks.test(merge.dat$`M+0`[merge.dat$EGFR=='Mut'],merge.dat$`M+0`[merge.dat$EGFR=='WT'])
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='WT']), col='#ccc6d6', main='Adjacent Lung',yaxt='n')
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']), pch=19, cex=1.5,
     xlim=c(50,100),add=T,main=NA,col="#df65b0",col.01line = "#df65b0")
lines(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']),col="#df65b0",pch=19,cex=0.6,
      col.points=brewer.pal(name = 'Set1',n = 6)[match(merge.dat$Pt.Code[merge.dat$EGFR=='Mut'],
                                                       unique(merge.dat$Pt.Code[merge.dat$EGFR=='Mut']))])
legend(68,1,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
title(ylab = "Cumulative Probability",
      outer = TRUE, line = 2.5,cex.lab=1.5)
title(xlab = "citrate m+0                 ",
      outer = TRUE, line = 2.5,cex.lab=1.5)

plot.new()
legend('center',pch=19,col=c('#df65b0','#ccc6d6'),legend = c('Mutant','WT'),bty='n',title = 'EGFR',xpd = T)

4H

cell_line<-sapply(lkb1$X,FUN = function(x) unlist(strsplit(x,' ',fixed = T))[1])
treatment<-sapply(lkb1$X,FUN = function(x) unlist(strsplit(x,' ',fixed = T))[2])
pv<-signif(sapply(unique(cell_line),FUN=function(cell) t.test(lkb1$M.0[cell_line==cell & treatment=='CTL'],lkb1$M.0[cell_line==cell & treatment=='LKB'])$p.value),2)
lkb1<-rbind(apply(lkb1[,-1],2,FUN=function(x) by(x,INDICES = lkb1[,1],mean)),
            apply(lkb1[,-1],2,FUN=function(x) by(x,INDICES = lkb1[,1],std.error)))
lkb1<-data.frame(X=rownames(lkb1),lkb1,X.1=rep(c('mean','SEM'),each=6))
lkb1.mean<-lkb1[lkb1$X.1=='mean',2:8]
lkb1.sem<-lkb1[lkb1$X.1=='SEM',2:8]
rownames(lkb1.mean)<-rownames(lkb1.sem)<-lkb1$X[1:(nrow(lkb1)/2)]
type<-rep('',nrow(lkb1.mean))
type[grep('CTL',rownames(lkb1.mean))]<-'CTL'
type[grep('OE',rownames(lkb1.mean))]<-'LKB1 OE'
lkb1.mean$Type<-factor(type)
lkb1.mean$`CellLine`<-factor(sapply(rownames(lkb1.mean),FUN=function(x) unlist(strsplit(x,split=' ',fixed = T))[1]))
tmp<-cbind(melt(lkb1.mean),melt(lkb1.sem))[,c(1:4,6)]
colnames(tmp)[4:5]<-c('mean','sem')
dat_text <- data.frame(
  x = c(1,1,1),
  y = c(57,57,57),
  label = pv,
  CellLine=levels(tmp$CellLine)
)
tmp$variable<-gsub('M.','m+',tmp$variable,fixed = T)
ggplot(tmp,aes(x=variable,y=mean,fill=Type))+ylab('% Enrichment')+
  geom_bar(stat="identity", position=position_dodge(), colour="black")+facet_wrap(~CellLine,ncol=1) +
  geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem), width=.2,
                position=position_dodge(.9))+theme_light()+ 
  scale_fill_manual("Cell Line", values = c("CTL" = "black", "LKB1 OE" = "white"))+ 
  annotate("rect", xmin = 0.75, xmax = 1.25, ymin = 50, ymax =50, alpha=1,colour = "black")+
  xlab(expression('citrate from [U-'^{13}*'C] glucose labeling'))+
  geom_text(data = dat_text,mapping = aes(x = x, y = y, label = label),inherit.aes = FALSE,size=2.5)

Fig 5

5C

AUC.erlotinib<-t(AUC[which(AUC$Common_chemical_name=='Erlotinib'),-(1:7)])
rownames(AUC.erlotinib)<-gsub('Calu','Calu-',rownames(AUC.erlotinib))
tmp<-merge(dat[,c('CitG6m0','CitQ6m5')],
           S4.data[,c('EGFR-pY1173','E-cadherin','Beta-catenin')],
           by.x=0,by.y = 0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
tmp<-merge(tmp,AUC.erlotinib,
           by.x=0,by.y = 0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
colnames(tmp)[6]<-'Erlotinib AUC'
tmp.5C<-tmp
tmp.col<-rep('black',nrow(tmp))
tmp.col[rownames(tmp)%in%c('H1869','H1650','HCC2935','HCC1897','HCC4019')]<-'purple'
tmp.col[rownames(tmp)%in%c('H1299','H1355','DFCI032','H1975','H23')]<-'green'
tmp$col<-factor(tmp.col)
colnames(tmp)<-gsub(' ','_',colnames(tmp),fixed = T)
colnames(tmp)<-gsub('-','_',colnames(tmp),fixed = T)
g<-ggpairs(tmp,
           upper=list(continuous=cor_fun),
           lower=list(mapping = aes(color=col,alpha=0.8),size=0.05),
           columns = 1:6)+theme_bw()+theme(panel.grid = element_blank())
colors<-c('darkgray','#4d9221','#c51b7d')
for(i in 2:6){
  for(j in 1:(i-1)){
    g[i,j]<-g[i,j] + scale_color_manual(values=colors)
  }
}
g

5E-figure

col.grad<-function(cols=c('blue','gray','red'),x){
  y=scales::alpha(colorRampPalette(cols)(length(x)),alpha = 1)[rank(x)]
  y[is.na(x)]<-'gray'
  return(y)
}

color.bar <- function(lut, min, max=-min, ticks=seq(min, max, len=2), title='') {
  scale = (length(lut)-1)/(max-min)
  plot(c(min,max), c(0,1), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='', main=title,mar=c(0,3,2,3),ylim=c(0,4))
  axis(3, ticks, las=1,labels = c('low','high'))
  for (i in 1:(length(lut)-1)) {
    x = (i-1)/scale + min
    rect(x,3,x+1/scale,4, col=lut[i], border=NA)
  }
}
tmp<-tmp.5C
colnames(tmp)[6]<-'AUC'
tmp$Name<-rownames(tmp)
tmp.h<-tmp[rownames(tmp)%in%c('H1869','H1650','HCC2935','HCC1897','HCC4019'),]
tmp.l<-tmp[rownames(tmp)%in%c('H1299','H1355','DFCI032','H1975','H23'),]
mut.shape=ifelse(mutations[match(rownames(tmp),rownames(mutations)),'EGFR']=="",19,2)
mut.shape[!is.na(str_extract(pattern = "p.[A-Z]7[0-9]+_[A-Z]7[0-9]+del",
                             string = mutations[match(rownames(tmp),rownames(mutations)),'EGFR']))]<-17
AUC.size<-rep(3,nrow(tmp))
AUC.size[is.na(AUC.erlotinib[match(rownames(tmp),rownames(AUC.erlotinib))])]<-1.5
ggplot(tmp, aes(x= CitG6m0, y = CitQ6m5, label=Name)) + 
  geom_point(
    colour=col.grad(cols=rev(brewer.pal(7,'PuOr')),
                    x=tmp$AUC),
    size=AUC.size,
    shape=mut.shape
  ) +
  geom_text_repel(data          = tmp.h,
                  colour='#c51b7d',
                  nudge_x       = -tmp.h$CitG6m0,
                  segment.size  = 0.2,
                  segment.color = "grey50") +
  geom_text_repel(data          = tmp.l,
                  colour="#4d9221",
                  nudge_x       = 70-tmp.l$CitG6m0,
                  segment.size  = 0.2,
                  segment.color = "grey50") +
  theme_classic(base_size = 15)

5E-legend

color.bar(colorRampPalette(rev(brewer.pal(7,'PuOr')))(100), -1)
legend('bottomleft',pch=c(19,17,2),legend = c('WT','exon19.del','Other'),title = 'EGFR status',bty = 'n')
legend('bottomright',pch=19,pt.cex = c(1,0.5),legend = c("available","missing"),col=c('black','grey'),title = 'Erlotinib AUC',bty = 'n')

5D

ge.emt<-ge[intersect(rownames(dat)[complete.cases(dat$CitG6m0)],rownames(ge)),intersect(colnames(ge),emt.signature)]
ge.emt[]<-apply(ge.emt,2,scale)
col.mat<-apply(dat[rownames(ge.emt),c('CitG6m0','CitQ6m5')],2,
               FUN=function(x) col.grad(cols=rev(brewer.pal(7,'PiYG')),x=x))
heatmap3(ge.emt,
         RowSideColors = col.mat,
         dist=dist,method = 'ward.D2',labCol = NA,labRow = NA)

5I-preparation and model fitting

tmp<-tmp.5C
colnames(tmp)[6]<-'AUC'
emt.pca<-prcomp(ge.emt)
mc<-Mclust(emt.pca$x[,1],G=2)$classification
tmp$`EMT-Signature`<-mc[match(rownames(tmp),rownames(ge.emt))]
tmp<-tmp[complete.cases(tmp),]
fit<-lm(AUC~.,tmp)
step<-stepAIC(fit, direction="both")
## Start:  AIC=446.85
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `Beta-catenin` + 
##     `EMT-Signature`
## 
##                   Df Sum of Sq    RSS    AIC
## - `Beta-catenin`   1      1212 215542 445.14
## - CitQ6m5          1      1911 216240 445.31
## - `EMT-Signature`  1      2807 217137 445.53
## <none>                         214329 446.85
## - `E-cadherin`     1     10665 224995 447.37
## - CitG6m0          1     16888 231218 448.79
## - `EGFR-pY1173`    1     50876 265205 455.92
## 
## Step:  AIC=445.14
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `EMT-Signature`
## 
##                   Df Sum of Sq    RSS    AIC
## - CitQ6m5          1      1639 217181 443.54
## - `EMT-Signature`  1      2293 217835 443.69
## <none>                         215542 445.14
## - `E-cadherin`     1      9883 225425 445.47
## - CitG6m0          1     15682 231224 446.79
## + `Beta-catenin`   1      1212 214329 446.85
## - `EGFR-pY1173`    1     50231 265773 454.04
## 
## Step:  AIC=443.54
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin` + `EMT-Signature`
## 
##                   Df Sum of Sq    RSS    AIC
## - `EMT-Signature`  1      2890 220071 442.22
## <none>                         217181 443.54
## - `E-cadherin`     1     11954 229135 444.32
## + CitQ6m5          1      1639 215542 445.14
## + `Beta-catenin`   1       941 216240 445.31
## - CitG6m0          1     21736 238917 446.50
## - `EGFR-pY1173`    1     48735 265915 452.06
## 
## Step:  AIC=442.22
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin`
## 
##                   Df Sum of Sq    RSS    AIC
## <none>                         220071 442.22
## - `E-cadherin`     1     11723 231794 442.92
## + `EMT-Signature`  1      2890 217181 443.54
## + CitQ6m5          1      2236 217835 443.69
## + `Beta-catenin`   1       424 219647 444.12
## - CitG6m0          1     19444 239515 444.63
## - `EGFR-pY1173`    1     50407 270478 450.95

5I-stepwise feature selection result

step$anova
## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `Beta-catenin` + 
##     `EMT-Signature`
## 
## Final Model:
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin`
## 
## 
##                Step Df Deviance Resid. Df Resid. Dev      AIC
## 1                                      45   214329.3 446.8493
## 2  - `Beta-catenin`  1 1212.298        46   215541.6 445.1426
## 3         - CitQ6m5  1 1639.044        47   217180.7 443.5365
## 4 - `EMT-Signature`  1 2890.472        48   220071.1 442.2240

5I-table

fit1<-lm(AUC~`CitG6m0`+`EGFR-pY1173`+`E-cadherin`,data = tmp)
fit2<-lm(AUC~`CitG6m0`+`EMT-Signature`+`EGFR-pY1173`+`E-cadherin`,data=tmp)
tab_model(fit1,fit2,
          dv.labels = c('Model 1 (AIC selected)','Model 2 (added EMT-Signature)'))
  Model 1 (AIC selected) Model 2 (added EMT-Signature)
Predictors Estimates CI p Estimates CI p
(Intercept) 553.43 503.88 – 602.98 <0.001 525.98 441.71 – 610.25 <0.001
Cit G 6 m 0 -1.89 -3.69 – -0.09 0.045 -2.04 -3.89 – -0.20 0.035
EGFR-pY1173 -28.06 -44.64 – -11.47 0.002 -27.64 -44.32 – -10.96 0.002
E-cadherin -11.22 -24.97 – 2.53 0.116 -18.26 -40.52 – 3.99 0.114
EMT-Signature 28.71 -42.43 – 99.84 0.433
Observations 52 52
R2 / adjusted R2 0.457 / 0.423 0.464 / 0.418

5J

yy<-data.frame(x=fit1$fitted.values,y=tmp$AUC)
ggplot(yy,aes(x=x,y=y))+geom_point(colour='navy',alpha=0.4,size=3)+
  geom_smooth(method='lm',se=F,colour=scales::alpha('black',0.3))+
  xlab('fitted values')+ylab('measured values')+ggtitle('Erlotinib AUC')+theme_classic()

Fig 6

6A

c<-rcorr(as.matrix(dat[,c('SerG6m3','GlyG6m2')]))
ggplot(dat[,c('SerG6m3','GlyG6m2')],mapping = aes(x=SerG6m3,y=GlyG6m2))+
         geom_point(colour='#6e016b',alpha=0.5,size=2.5)+
         theme_bw()+
         ggtitle(paste0("Serine Glycine interconversion\nr = ",round(c$r[1,2],2),', pv = ',signif(c$P[1,2],2)))

6B

drug<-S4.data[,grep(")$",colnames(S4.data))]
rownames(drug)<-S4.data$Cell.Line
drug<-as.matrix(drug[match(rownames(dat),rownames(drug)),]);class(drug)<-'numeric'
pv<-apply(drug,2,FUN=function(x) cor.test(dat$SerG6m3,x,method = 'pearson')$p.value)
r<-apply(drug,2,FUN=function(x) cor.test(dat$SerG6m3,x,method = 'pearson')$estimate)
n<-apply(drug,2,FUN=function(x) length(which(!is.na(x) & !is.na(dat$SerG6m3))))
res<-data.frame(pv,r,n)
rownames(res)<-sapply(rownames(res),FUN=function(x) unlist(strsplit(x,split = '.',fixed = T))[1])
res<-res[order(res$pv,decreasing = F),]
res$pv<-(-log10(res$pv))
res$drug<-factor(rownames(res),levels = rev(rownames(res)))
res$significance<-rep('insig',nrow(res))
res$significance[1:countSignificant(Bum(pv),0.05)]<-'sig'
res$significance<-factor(res$significance,levels = c('insig','sig'))
ggplot(res,aes(y=pv,fill=significance,x=drug))+geom_bar(stat="identity")+coord_flip() +scale_fill_brewer(palette="PuRd")+
  geom_hline(yintercept = (-log10(0.05)), linetype="dotted")+theme_classic()+ylab('-log10(p-value)')+xlab("")+ggtitle('Correlation with SerG6m3')

6D

tmp<-data.frame(Pemetrexed=drug[,"Pemetrexed.(uM)"],SerG6m3=dat$SerG6m3)
rownames(tmp)<-rownames(dat)
tmp<-tmp[complete.cases(tmp),]
mc<-Mclust(tmp$Pemetrexed,G=2)
tmp$sensitivity<-c('sensitive','resistant')[mc$classification]
tmp$sensitivity<-factor(tmp$sensitivity,levels=c('sensitive','resistant'))
tmp$cutoff<-c('low','high')[as.numeric(tmp$SerG6m3>16)+1]
h.cells<-c('PC-9','H2170','H596','H460','H1299')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195')
tmp$name<-rownames(tmp)
tmp.h<-tmp[rownames(tmp)%in%h.cells,]
tmp.l<-tmp[rownames(tmp)%in%l.cells,]
h.col<-brewer.pal(n = 3,'Set1')[1]
l.col<-brewer.pal(n = 3,'Set1')[2]
ggplot(tmp,aes(x=SerG6m3,y=Pemetrexed,colour=sensitivity,shape=cutoff, label=name))+geom_point()+
  geom_text_repel(data          = tmp.h,
                  colour=h.col,
                  nudge_y       = 1-tmp.h$Pemetrexed,
                  segment.size  = 0.2,
                  segment.color = "grey50") +
  geom_text_repel(data          = tmp.l,
                  colour=l.col,
                  nudge_x       = 70-tmp.l$SerG6m3,
                  segment.size  = 0.2,
                  segment.color = "grey50") +
  theme_bw()

6E

mat.mean<-function(x) apply(x,2,FUN=function(y) mean(y,na.rm = T))
pem.mean<-do.call(rbind,by(data = pem[,-1],INDICES = pem[,1],FUN = mat.mean))
pem.auc<-apply(pem.mean,2,FUN=function(x) trapz(log10(as.numeric(rownames(pem.mean))),x))
pem.plot<-data.frame(cell=names(pem.auc),auc=pem.auc)
h.cells<-c('PC-9','H2170','H596','H460','H1299')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195')
pem.plot$SerG6m3<-rep(NA,nrow(pem.plot))
pem.plot$SerG6m3[as.character(pem.plot$cell)%in%h.cells]<-'high'
pem.plot$SerG6m3[as.character(pem.plot$cell)%in%l.cells]<-'low'
pem.plot$SerG6m3<-factor(pem.plot$SerG6m3,levels = c('high','low'))
pem.plot$cell<-factor(pem.plot$cell,levels = c(h.cells,l.cells))
pem.plot$mean<-unlist(by(data = pem.plot$auc,INDICES = pem.plot$SerG6m3,mean))[match(pem.plot$SerG6m3,c('high','low'))]
pem.plot$sd<-unlist(by(data = pem.plot$auc,INDICES = pem.plot$SerG6m3,sd))[match(pem.plot$SerG6m3,c('high','low'))]
set.seed(621)
ggplot(pem.plot,aes(x=SerG6m3,y=auc,colour=SerG6m3))+ geom_jitter(width = 0.2)+scale_color_brewer(palette="Set1")+
  ggtitle(paste0('pv = ',signif(t.test(auc~SerG6m3,pem.plot)$p.value,1)))+ylab('AUC')+
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd, group = SerG6m3),width = 0.2)+theme_bw()

6F-H

g.list<-list()
for(sheet in sheets(wb.misty)){
  tmp<-readWorkbook(wb.misty,sheet,rowNames = T)
  tmp<-data.frame(days=rownames(tmp),Cell=sheet,melt(tmp))
  tmp$variable<-sapply(as.character(tmp$variable),FUN=function(x) substr(x,1,nchar(x)-2))
  tmp$variable<-factor(tmp$variable,levels = c('Control','Pem.500', 'Pem.750','Pem.1000'))
  tmp<-tmp[complete.cases(tmp),]
  tmp$days<-as.character(tmp$days)
  colnames(tmp)[colnames(tmp)=='variable']<-'Treatment'
  g.list[[sheet]]<-ggpar(ggline(tmp[tmp$Cell==sheet,], x = "days", y = "value", add = "mean_se",color = 'Treatment',
                                palette = "Dark2")+ylab('tumor volume')+
                           stat_compare_means(aes(group = Treatment), label = "p.signif", method = 'aov',
                                              label.y = rep(max(na.omit(unlist(by(data = tmp$value,INDICES = list(tmp$Treatment,tmp$days),mean))))*1.2,
                                                            length(unique(tmp$days[tmp$Cell==sheet]))))+
                           ggtitle(sheet),legend='none')
}
marrangeGrob(g.list, nrow=1, ncol=3,top=NULL)

6F-H-legend

legend <- cowplot::get_legend(ggpar(ggline(tmp[tmp$Cell==sheet,], x = "days", y = "value", add = "mean_se",color = 'Treatment',
                              palette = "Dark2")+ylab('tumor volume')+
                         stat_compare_means(aes(group = Treatment), label = "p.signif", method = 'aov',
                                            label.y = rep(max(na.omit(unlist(by(data = tmp$value,INDICES = list(tmp$Treatment,tmp$days),mean))))*1.2,
                                                          length(unique(tmp$days[tmp$Cell==sheet]))))+
                         ggtitle(sheet),legend='right'))
grid.newpage();grid.draw(legend)

Fig S1

S1A

doubling.time=(48/log2(dat$`Day3/Day1`)+96/log2(dat$`Day5/Day1`))/2
nutrients<-dat[,c('Lac','Glc','Lac/Glc','Glu','Gln','Glu/Gln')]
nutrients.norm<-nutrients
nutrients.norm[]<-apply(nutrients,2,FUN=function(x) x/((2^(7/doubling.time)+1)/2))
nutrients.norm$`Glu/Gln`<-nutrients.norm$Glu/nutrients.norm$Gln
nutrients.norm$`Lac/Glc`<-nutrients.norm$Lac/nutrients.norm$Glc
tmp<-data.frame(melt(nutrients[,c(1,2,4,5)]),corrected=melt(nutrients.norm[,c(1,2,4,5)])[,2])
colnames(tmp)[1:2]<-c('feature','uncorrected')
tmp$cor<-rep(unlist(lapply(names(nutrients)[c(1,2,4,5)],FUN=function(x) round(cor.test(nutrients[,x],nutrients.norm[,x])$estimate,3))),each=nrow(nutrients))
tmp$feature_name=factor(paste(tmp$feature,tmp$cor,sep="\nr = "),ordered = F)
tmp$feature_name=factor(tmp$feature_name,levels(tmp$feature_name)[c(1,4,2,3)])
ggplot(tmp,aes(x=uncorrected,y=corrected))+geom_point(alpha=0.3)+theme_classic()+facet_wrap(~feature_name,scale='free',ncol = 1)+
  xlab('normalized to final protein content')+ylab('normalized to estimated average protein content in 7h growth')

S1B

replicate.sd<-list()
cell.sd<-list()
new.names<-c('CitG6','CitG24','CitQ6','CitQ24','SerG6','SerG24',
             'GlyG6','GlyG24','FumG6','FumG24','FumQ6','FumQ24',
             'MalG6','MalG24','MalQ6','MalQ24','LacG6','LacG24','LacQ6','LacQ24')
result.list<-list()
sd.mat.list<-list()
for(i in 2:(length(new.names)+1)){
  sheet<-sheets(wb)[i]
  tmp<-read.xlsx(wb, sheet = sheet,rowNames = T,colNames = T)
  sd.mat<-as.matrix(tmp[-1,grep('Std',colnames(tmp),ignore.case = T):(grep('Rep',colnames(tmp),ignore.case = T)-1)])
  rownames(sd.mat)<-tmp[-1,1]
  sd.mat<-sd.mat[complete.cases(sd.mat),]
  class(sd.mat)<-'numeric'
  btw.cell.sd<-apply(as.matrix(tmp[-1,grep('ave',colnames(tmp),ignore.case = T):(grep('Std',colnames(tmp),ignore.case = T)-1)]),2,FUN=function(x) sd(x,na.rm = T))
  
  result.list[[length(result.list)+1]]<-data.frame(feature=rep(sheet,length(btw.cell.sd)),btw.cell.sd,btw.rep.sd.mean=apply(sd.mat,2,mean),btw.rep.sd.sd=apply(sd.mat,2,sd))
  sd.mat.list[[length(sd.mat.list)+1]]<-sd.mat
  names(sd.mat.list)[length(sd.mat.list)]<-new.names[i-1]
}
tmp<-do.call(rbind,result.list)
tmp$metabolite<-substr(tmp$feature,1,3)
tmp$`labeling time`<-factor(paste(substr(tmp$feature,5,8),'hr'),levels=c('6 hr','24 hr'))
tmp$tracer<-c('glucose','glutamine')[match(substr(tmp$feature,4,4),c('G','Q'))]

ggplot(tmp,aes(x=btw.cell.sd,y=btw.rep.sd.mean,colour=metabolite,shape=tracer,alpha=`labeling time`))+
  geom_abline(aes(intercept=0,slope=1),colour='lightgray')+
  scale_alpha_discrete(range=c(0.3,0.8))+
  geom_point()+coord_flip()+
  geom_pointrange(aes(ymin=btw.rep.sd.mean, ymax=btw.rep.sd.mean+btw.rep.sd.sd))+
  theme_classic()+ylim(0,22)+xlim(0,22)

S1C was made by Pei-Hsuan Chen

S1D

diversity.plot<-function(sheet,cell_lines,tracer){
  tmp<-read.xlsx(wb, sheet = sheet,rowNames = T,colNames = T,startRow = 2)
  tmp$Cell.line<-toupper(tmp$Cell.line)
  tmp<-tmp[tmp$Cell.line%in%cell_lines,]
  rownames(tmp)<-tmp$Cell.line;tmp<-tmp[,-c(1:ifelse(sheet=='CitG6',16,15))]
  tmp.mean<-apply(t(tmp),2,FUN=function(x){
    apply(matrix(byrow = F,data = x,nrow = 7,ncol = length(x/7)),1,FUN = function(y) mean(na.omit(y)))
  })
  tmp.sem<-apply(t(tmp),2,FUN=function(x){
    apply(matrix(byrow = F,data = x,nrow = 7,ncol = length(x/7)),1,FUN = function(y) std.error(na.omit(y)))
  })
  rownames(tmp.sem)<-rownames(tmp.mean)<-paste('m+',0:6,sep = '')
  tmp.m<-setNames(cbind(melt(tmp.mean),melt(tmp.sem)[,3]),c('isotopologue','cell line','mean','sem'))
  if(tracer=='glucose'){
    title<-expression('citrate labeling by [U-'^{13}*'C] glucose labeling')
  }else{
    title<-expression('citrate labeling by [U-'^{13}*'C] glutamine labeling')
  }
  ggplot(tmp.m,aes(x=isotopologue,y=mean))+
    geom_bar(stat="identity",position=position_dodge())+
    geom_errorbar(aes(ymin=mean-sem,ymax=mean+sem),width=.2,position=position_dodge(.9))+
    ylab('% enrichment')+
    ggtitle(title)+
    theme_bw()+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))+
    facet_grid(.~`cell line`)
}
diversity.plot(sheet='CitG6',cell_lines = c('HCC44','H1869','H1395'),tracer='glucose')

S1E

diversity.plot(sheet='CitQ6',cell_lines = c('H1975','H441','H1395'),tracer='glutamine')

Fig S2

S2A

Isotopologue.Distribution(Key = 'Fum',Name='Fumarate')

S2B

Isotopologue.Distribution(Key = 'Mal',Name='Malate')

S2C

Isotopologue.Distribution(Key = 'Lac',Name='Lactate')

S2D

Isotopologue.Distribution(Key = 'Ser',Name='Serine')

S2E

Isotopologue.Distribution(Key = 'Gly',Name='Glycine')

S2F

dat.tracing<-dat[,grep('m',colnames(dat))]
dat.tracing.var<-apply(dat.tracing,2,FUN=function(x) sd(x,na.rm = T))
var.df<-data.frame(variance=dat.tracing.var,
                   metabolite=factor(substr(names(dat.tracing.var),1,3),levels = c('Cit','Fum','Mal','Lac','Ser','Gly')),
                   timepoint=factor(ifelse(grepl('24',names(dat.tracing.var)),'24h','6h'),levels=c('6h','24h')),
                   tracer=factor(ifelse(grepl('Q',names(dat.tracing.var)),'Glutamine','Glucose'),levels = c('Glucose','Glutamine')),
                   isotopologue=as.factor(paste('m+',sapply(names(dat.tracing.var),FUN=function(x) rev(unlist(strsplit(x,split = '')))[1]))))
ggplot(data = var.df,
       aes(x=isotopologue,y=variance,fill=timepoint))+
  geom_bar(stat="identity")+
  facet_wrap(metabolite~tracer+timepoint,ncol = 4, scales = "free_x")+theme_light()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Fig S3

dat.6h<-dat[,grep('24',colnames(dat),invert = T)]
colnames(dat.6h)<-gsub("G6","G",colnames(dat.6h))
colnames(dat.6h)<-gsub("Q6","Q",colnames(dat.6h))
corrplot(corr = cor(dat.6h,use = 'complete.obs'),type= 'upper' ,diag = F,col=col2,tl.col = 'black',tl.cex = 0.8)

Fig S4

S4B-D

tmp<-merge(ge[,'PC',drop=F],dat[,c('CitG6m4','MalG6m3')],by.x=0,by.y=0,all=T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
tmp<-tmp[complete.cases(tmp),]
fit<-lm(MalG6m3~CitG6m4,tmp)
tmp$res<-lm(MalG6m3~CitG6m4,tmp)$residuals
par(mfrow=c(2,2))
plot.new()
plot(tmp[,c('CitG6m4','MalG6m3')],type='n')
abline(lm(MalG6m3~CitG6m4,tmp))
for(cell in names(fit$residuals)){
  segments(x0 = tmp$CitG6m4[rownames(tmp)==cell],x1=tmp$CitG6m4[rownames(tmp)==cell],
           y0 = tmp$MalG6m3[rownames(tmp)==cell],y1 = fit$fitted.values[cell],col='#f68712')
}
points(tmp[,c('CitG6m4','MalG6m3')],pch=19)
legend('topleft',col='#f68712',legend = 'residual',lty=1)
ct<-cor.test(tmp$MalG6m3,tmp$PC)
plot(tmp$MalG6m3,tmp$PC,pch=19,xlab='MalG6m3',ylab='PC gene expression',
     main=paste('r = ',round(ct$estimate,2),'; pv = ',signif(ct$p.value,2),sep = ''))
ct<-cor.test(tmp$res,tmp$PC)
plot(tmp$res,tmp$PC,pch=19,xlab='residuals',ylab='PC gene expression',
     main=paste('r = ',round(ct$estimate,2),'; pv = ',signif(ct$p.value,2),sep = ''))

Fig S5

S5C

h<-hclust(dist(ge.emt),method = 'ward.D2')
hc<-cutree(h,2)
# the line commented below was used to check which cluster was the epithelial cluster
# boxplot(ge.emt[,'CDH1']~hc)
E<-rownames(ge.emt)[hc==2]
M<-rownames(ge.emt)[hc==1]
metab.p<-unlist(lapply(c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln'),
                       FUN=function(x) t.test(na.omit(dat[rownames(dat)%in%E,x]),na.omit(dat[rownames(dat)%in%M,x]))$p.value))
names(metab.p)<-c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')
tmp<-data.frame(dat[match(c(E,M),rownames(dat)),c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')],
                emt=rep(c('epithelial','mesenchymal'),c(length(E),length(M))),check.names = F)
dat_text <- data.frame(
  x = rep(1.85,4),
  y = 0.9*apply(tmp[,1:4],2,FUN=function(x) max(x,na.rm = T)),
  label = paste("pv =",signif(metab.p,2)),
  variable=c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')
)
tmp<-melt(tmp)
ggplot(tmp,aes(x=emt,y=value))+geom_boxplot()+
  geom_text(data = dat_text,mapping = aes(x = x, y = y, label = label),inherit.aes = FALSE,size=4)+
  facet_wrap(~variable,scale='free_y')+xlab("")+theme_bw()

S5D

tmp<-merge(CDH1.data,dat[,'CitG6m0',drop=F],by.x=0,by.y=0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
tmp<-tmp[complete.cases(tmp),]
ggpairs(tmp[complete.cases(tmp),],
        upper=list(continuous=cor_fun),
        lower=list(mapping = aes(alpha=0.8),size=0.05))+
  theme_bw()+theme(panel.grid = element_blank())

S6A

common<-intersect(rownames(dat)[!is.na(dat$SerG6m3)],rownames(ge))
dat.Ser<-dat[common,]
ge.Ser<-ge[common,c('PHGDH','PSAT1','PSPH','CBS','SHMT1','SHMT2')]
ge.Ser<-ge.Ser[order(dat.Ser$SerG6m3,decreasing = T),]
ge.Ser[]<-apply(ge.Ser,2,rank)

annot_df<-dat.Ser[order(dat.Ser$SerG6m3,decreasing = T),'SerG6m3',drop=F]
annot_df[,1]<-rank(annot_df[,1])
mf_range<-range(annot_df[,1])
col<-circlize::colorRamp2(seq(from=mf_range[1],to=mf_range[2],length.out = 3), 
                          c(l.col, "#EEEEEE", h.col))
col<-list(col)
names(col)<-'SerG6m3'
# Create the heatmap annotation
ha <- HeatmapAnnotation(annot_df, col = col, show_legend = F)
colnames(ge.Ser)<-paste(colnames(ge.Ser),
                           "\n",'rho=',round(apply(ge.Ser,2,FUN=function(x) cor.test(x,sort(dat.Ser$SerG6m3,decreasing = T),method = 'spearman')$estimate),1),
                           ', p=',signif(apply(ge.Ser,2,FUN=function(x) cor.test(x,sort(dat.Ser$SerG6m3,decreasing = T),method = 'spearman')$p.value),1),sep = '')

ht<-Heatmap(t(ge.Ser), 
            show_row_names = T,show_column_names = F,
            cluster_rows = F,cluster_columns = F,
            show_heatmap_legend = F,
            top_annotation = ha)
lgd = list(Legend(at = c("low","","medium","","high"), title = "gene expr", type = "grid", 
                  legend_gp = gpar(fill = c("#0000FF","#7777F6","#EEEEEE","#F67777","#FF0000"))),
           Legend(at = c("low",'medium',"high"), title = "SerG6m3", type = "grid", 
                  legend_gp = gpar(fill = c(l.col, "#EEEEEE", h.col))))

draw(ht, column_title = 'Serine metabolism genes',annotation_legend_list = lgd)
an<-'SerG6m3'
decorate_annotation(an, {
  # annotation names on the left
  grid.text(an, unit(1, "npc") + unit(2, "mm"), 0.5, default.units = "npc", just = "left")
})

S6B

tmp.cor<-rcorr(as.matrix(dat[,c('SerG6m3','Day3/Day1')]))
ggplot(data = dat[,c('SerG6m3','Day3/Day1')],aes(x=SerG6m3,y=`Day3/Day1`))+geom_point(color=scales::alpha('#8856a7',0.7))+
  ggtitle(paste0('r = ',round(tmp.cor$r[1,2],2),', pv = ',signif(tmp.cor$P[1,2],2)))+theme_bw()

S6C-J

h.cells<-c('PC-9','H2170','H596','H460','H1299','H1155')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195','H2250')
hl.plot<-function(tmp,feature,hide=T){
  tmp.plt<-data.frame(cell=names(tmp),value=tmp)
  tmp.plt$SerG6m3<-rep(NA,nrow(tmp.plt))
  tmp.plt$SerG6m3[tmp.plt$cell%in%h.cells]<-'high'
  tmp.plt$SerG6m3[tmp.plt$cell%in%l.cells]<-'low'
  tmp.plt$SerG6m3<-factor(tmp.plt$SerG6m3,levels = c('high','low'))
  tmp.plt$mean<-unlist(by(data = tmp.plt$value,INDICES = tmp.plt$SerG6m3,mean))[match(tmp.plt$SerG6m3,c('high','low'))]
  tmp.plt$sd<-unlist(by(data = tmp.plt$value,INDICES = tmp.plt$SerG6m3,sd))[match(tmp.plt$SerG6m3,c('high','low'))]
  tmp.plt<-tmp.plt[complete.cases(tmp.plt),]
  pv<-t.test(value~SerG6m3,tmp.plt)$p.value
  g<-ggplot(tmp.plt,aes(x=SerG6m3,y=value,colour=SerG6m3))+ geom_jitter(width = 0.2)+scale_color_brewer(palette="Set1")+
    ggtitle(paste('pv =',signif(pv,2)))+ylab(feature)+
    geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd, group = SerG6m3),width = 0.2)+theme_bw()
  if(hide) g+theme(legend.position = "none")
  else g
}


g.list<-list()
g.list[[1]]<-hl.plot(apply(brdu,2,mean),feature='Brdu Staining')
for(i in 1:4){
  g.list[[i+1]]<-hl.plot(apply(as.matrix(nova[[i]]),2,mean),feature=names(nova)[i])
}
g.list[[6]]<-hl.plot(apply(ser.pool,2,mean),feature='Serine Pool Size')
g.list[[7]]<-hl.plot(apply(as.matrix(pt.all$Ser),2,mean),feature='Intracellular Serine')
g.list[[8]]<-hl.plot(apply(as.matrix(med.all$Ser),2,mean),feature='Extracellular Serine')
marrangeGrob(g.list, layout_matrix = matrix(1:9, byrow = T, nrow = 3),top=NULL)

S6C-J-legend

legend <- cowplot::get_legend(hl.plot(apply(brdu,2,mean),feature='Brdu Staining',hide = F))
grid.newpage();grid.draw(legend)

sessionInfo

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server 7.5 (Maipo)
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2        cowplot_0.9.4         ggpubr_0.2           
##  [4] magrittr_1.5          pracma_2.2.2          ClassComparison_3.1.6
##  [7] oompaBase_3.2.6       ComplexHeatmap_1.12.0 sjPlot_2.6.2         
## [10] mclust_5.4.2          plotrix_3.7-4         stringr_1.4.0        
## [13] heatmap3_1.1.1        ggrepel_0.8.0.9000    factoextra_1.0.5     
## [16] gridExtra_2.3         gtable_0.3.0          reshape2_1.4.3       
## [19] ppcor_1.1             MASS_7.3-51           GSVA_1.22.4          
## [22] RColorBrewer_1.1-2    corrplot_0.84         Hmisc_4.1-1          
## [25] Formula_1.2-3         survival_2.42-6       lattice_0.20-35      
## [28] openxlsx_4.1.0        GGally_1.4.0          ggplot2_3.1.0        
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.3      circlize_0.4.5       igraph_1.2.2        
##   [4] plyr_1.8.4           lazyeval_0.2.1       GSEABase_1.36.0     
##   [7] TMB_1.7.14           splines_3.3.2        TH.data_1.0-10      
##  [10] digest_0.6.18        htmltools_0.3.6      viridis_0.5.1       
##  [13] checkmate_1.8.5      memoise_1.1.0        cluster_2.0.7-1     
##  [16] fastcluster_1.1.25   annotate_1.52.1      modelr_0.1.3        
##  [19] sandwich_2.5-0       colorspace_1.3-2     blob_1.1.1          
##  [22] haven_2.0.0          xfun_0.7             dplyr_0.7.6         
##  [25] crayon_1.3.4         RCurl_1.95-4.11      graph_1.52.0        
##  [28] lme4_1.1-20          bindr_0.1.1          zoo_1.8-4           
##  [31] glue_1.3.0           emmeans_1.3.2        GetoptLong_0.1.7    
##  [34] sjstats_0.17.3       sjmisc_2.7.7         kernlab_0.9-27      
##  [37] shape_1.4.4          prabclus_2.2-7       BiocGenerics_0.20.0 
##  [40] DEoptimR_1.0-8       scales_0.5.0         mvtnorm_1.0-8       
##  [43] DBI_1.0.0            ggeffects_0.8.0      Rcpp_0.12.17        
##  [46] viridisLite_0.3.0    xtable_1.8-3         htmlTable_1.12      
##  [49] foreign_0.8-71       bit_1.1-14           stats4_3.3.2        
##  [52] prediction_0.3.6.2   htmlwidgets_1.3      fpc_2.1-11.1        
##  [55] acepack_1.4.1        modeltools_0.2-22    pkgconfig_2.0.2     
##  [58] reshape_0.8.8        XML_3.98-1.16        flexmix_2.3-14      
##  [61] nnet_7.3-12          labeling_0.3         tidyselect_0.2.4    
##  [64] rlang_0.2.1          AnnotationDbi_1.36.2 munsell_0.5.0       
##  [67] tools_3.3.2          generics_0.0.2       RSQLite_2.1.1       
##  [70] sjlabelled_1.0.16    broom_0.5.1          ggridges_0.5.1      
##  [73] evaluate_0.13        yaml_2.2.0           knitr_1.23          
##  [76] bit64_0.9-7          zip_1.0.0            robustbase_0.93-3   
##  [79] purrr_0.2.5          dendextend_1.9.0     coin_1.2-2          
##  [82] nlme_3.1-128         whisker_0.3-2        bayesplot_1.6.0     
##  [85] rstudioapi_0.9.0     tibble_1.4.2         stringi_1.2.4       
##  [88] forcats_0.3.0        trimcluster_0.1-2.1  Matrix_1.2-14       
##  [91] psych_1.8.12         nloptr_1.2.1         ggsci_2.9           
##  [94] stringdist_0.9.5.1   pillar_1.3.0         pwr_1.2-2           
##  [97] GlobalOptions_0.1.0  estimability_1.3     data.table_1.11.8   
## [100] bitops_1.0-6         R6_2.3.0             latticeExtra_0.6-28 
## [103] IRanges_2.8.2        codetools_0.2-15     assertthat_0.2.0    
## [106] rjson_0.2.20         withr_2.1.2          mnormt_1.5-5        
## [109] multcomp_1.4-8       S4Vectors_0.12.2     diptest_0.75-7      
## [112] parallel_3.3.2       hms_0.4.2            rpart_4.1-13        
## [115] tidyr_0.8.2          coda_0.19-2          glmmTMB_0.2.2.0     
## [118] class_7.3-14         minqa_1.2.4          rmarkdown_1.13      
## [121] snakecase_0.9.2      Biobase_2.34.0       base64enc_0.1-3